Accepted by ApJ: Jan 28, 2000 



Structure and Evolution of the Envelopes of Deeply Embedded Massive Young 

Stars 

Floris F. S. van der Tak & Ewine F. van Dishoeck 
, Sterrewacht, Postbus 9513, 2300 RA Leiden, The Netherlands 

o 



o 

(N 



O 

o 
o 



Neal J. Evans II 



. Department of Astronomy, University of Texas, Austin, TX 78712 

and 

CN ■ Geoffrey A. Blake 



Division of Geological and Planetary Sciences, California Institute of Technology, MS 150-21, 



in ■ Pasadena, CA 91125 

ABSTRACT 



Qh| The physical structure of the envelopes around a sample of fourteen massive young 

O ' stars is investigated using maps and spectra in submillimeter continuum and lines 

' of C^^O, CS, C^^S and H2CO. Nine of the sources are highly embedded luminous 

■ (10^ — 10^ Lq) young stellar objects which are bright near-infrared sources, but weak 

in radio continuum; the other objects are similar but not bright in the near-infrared, 
^ . and contain "hot core"-type objects and/or ultracompact H II regions. The data are 

used to constrain the temperature and density structure of the circumstellar envelopes 
on 10^ — 10^ AU scales, to investigate the relation between the different objects and to 
search for evolutionary effects. 

The total column densities and the temperature profiles are obtained by fitting 
self-consistent dust models to submillimeter photometry. The calculated temperatures 
range from 300 to 1000 K at ~ 10^ AU and from 10 to 30 K at ~ 10^ AU from the 
star. Visual extinctions are a few hundred to a few thousand magnitudes, assuming a 
grain opacity at A1300/im, of ~ 1 cm~^ g~^ of dust, as derived earlier for one of our 
sources. The mid-infrared data are consistent with a 30% decrease of the opacity at 
higher temperatures, caused by the evaporation of the ice mantles. 

The CS, C^'^S and II2CO data as well as the submillimeter dust emission maps 
indicate density gradients n cx r~°. Assuming a constant CS abundance throughout the 
envelope, values of a = 1.0 — 1.5 are found, significantly flatter than the a = 2.0 it 0.3 
generally found for low-mass objects. This flattening may indicate that in massive 
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young stellar objects, nonthermal pressure is more important for the support against 
gravitational collapse, while thermal pressure dominates for low-mass sources. Wc find 
a = 2 for two hot core-type sources, but regard this as an upper limit since in these 
objects, the CS abundance may be enhanced in the warm gas close to the star. 

The assumption of spherical symmetry is tested by modeling infrared absorption 
line data of ^^CO, CS emission line profiles and near-infrared continuum. There is a 
distinct, but small deviation from spherical symmetry: the data are consistent with a 
decrease of the optical depth by a factor of ~ 3 in the central < 10". The homogeneity 
of the envelopes is verified by the good agreement of the total masses in the power law 
models with the virial masses. 

Modeling of C^^O emission shows that 40 — 90% of the CO is frozen out onto the 
dust. The CO abundances show a clear correlation with temperature, as expected if 
the abundance is controlled by freeze-out and thermal desorption. The CS abundance 
is 3 X 10~^ on average, ranging from (4 — 8) x 10"^*^ in the cold source GL 7009S to 
(1 — 2) X 10^*^ in the two "hot core" -type sources. 

Dense outflowing gas is seen in the CS and H2CO line wings; the predominance of 
blueshifted emission suggests the presence of dense, optically thick material within 10" of 
the center. Interferometric continuum observations at A1300 — 3500/im show compact 
emission, probably from an 0'.'3 diameter, optically thick dust component, such as a 
dense shell or a disk. The emission is a factor of 10 — 100 stronger than expected for the 
envelopes seen in the single-dish data, so that this component may be opaque enough 
to explain the asymmetric CS and H2CO line profiles. 

The evolution of the sources is traced by the overall temperature (measured by 
the far- infrared color), which increases systematically with decreasing ratio of envelope 
mass to stellar mass. The observed anticorrelation of near-infrared and radio contin- 
uum emission suggests that the erosion of the envelope proceeds from the inside out. 
Conventional tracers of the evolution of low-mass objects do not change much over this 
narrow age range. 

Subject headings: ISM: dust, extinction, ISM: jets and outflows, ISM: molecules. Stars: 
circumstellar matter; Stars: Formation. 

1. Introduction 

The dynamical processes governing massive star formation are much less understood than is 
the case for low-mass stars (Churchwell 1999; Garay & Lizano 1999). The various observational 
"appearances" of low-mass star formation (T Tauri stars, FU Orionis stars, infrared Class I-III 
sources, molecular outflows, ...) have been linked into a single, if rough, general evolutionary 
scenario (Shu et al. 1993; Andre et al. 1993, 1999; Evans 1999). In contrast, no clear evolutionary 
sequence has been established for high-mass stars. Objects such as the BN object in Orion are 
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highly embedded and emit the bulk of their luminosity in the mid- and far-infrared. Ultracompact 

(UC) H II regions are small (~ 0.1 pc) sources of free-free emission at radio wavelengths (Wood & 
Churchwell 1989; Kurtz et al. 1994). "Hot core"-typc objects have bright molecular line emission 
at submillimeter wavelengths, which indicate temperatures of several 100 K and high abundances 
of saturated carbon-bearing molecules (e.g., Blake et al. 1987). They usually lack detectable radio 
emission, which might be due to "quenching" by material accreting at rates >10~^ Mq yr~^ 
(Walmsley 1995). The distinction between UC H II regions and hot cores is not always clear, 
and the two are often found located very close to each other (Cesaroni et al. 1999; Kurtz et al. 
1999). As a step towards understanding the evolution of massive YSOs, this paper presents models 
for a sample of fourteen such objects. 

The formation of high-mass stars is invisible at optical wavelengths because of the high opacity 
of the surrounding material, so that reliable age estimates for massive protostellar objects are diffi- 
cult to obtain. Such estimates for massive protostars have traditionally come from the morphology 
of the radio continuum emission (Colley 1980), but this method applies only to relatively evolved 
objects. More recently, the size to velocity ratio or dynamical time scale of molecular outflows has 
been used as a measure of age, but the derived age depends strongly on the adopted dynamical 
model (Cabrit et al. 1997). The chemical composition of the material is a potentially powerful 
clock, but difficult to calibrate (Charnley et al. 1992; van Dishoeck & Blake 1998). 

The density and temperature structure of the circumstellar material are also expected to 
change as the central star develops. For the density, power laws n(r) oc r~" are predicted by 
many theoretical models. Before collapse begins, low-mass cores may relax to a power law, with 
a = 2.0 for thermal support (Shu 1977), and a = 1.0 for turbulent support (Lizano & Shu 1989; 
Myers & Fuller 1992; McLaughlin & Pudritz 1996). Once collapse begins, the density distribution 
tends toward a = 1.5 (Larson 1969; Shu 1977). The temperature structure T(r) is determined by 
the luminosity, the dust opacities, and n(r). Since the response of T(r) to changes in luminosity 
is rapid, we use the observed luminosity to determine T(r). Besides giving information on the 
dynamics of star formation, a good model of the physical structure of objects is a prerequisite for 
determining their chemical composition, which then by itself can be used as a powerful additional 
evolutionary indicator. The physical conditions around newly-born massive stars also reveal some 
of the influence that such stars have on their surroundings, which is of interest for the effects of 
star formation on large ( ~ 1 pc) scales. 

The physical structure of massive YSO envelopes has been studied based on dust continuum 
observations by, e.g., Chini et al. (1986), Churchwell et al. (1990), Faison et al. (1998) and Cesaroni 
et al. (1998). Such studies are sensitive to the column density and the temperature of the envelope, 
but not to density itself. A major source of uncertainty for all dust models is the choice of optical 
properties of the grains, which limits the accuracy of derived envelope masses to a factor of 2 
(Henning 1997). 

Molecular rotational lines are direct probes of the H2 density and are also sensitive to temper- 
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ature. Zhou et al. (1991, 1994) and Carr et al. (1995) investigated the structure of massive star- 
forming regions using Hues of CS, H2CO and other molecules with large dipole moments. However, 
in these studies, the column density (or the cloud mass) was not independently constrained. Doing 
this requires knowledge of the molecular abundance, information which is generally unavailable 
except in the case of CO, which is relatively inert and much more abundant than most molecules. 
Still, part of the CO may be frozen out onto grain surfaces in the cold regions far from embedded 
stars. 

This paper uses observations of both the dust and the molecular gas to constrain the physical 
structure of the envelopes of a sample of massive young stars. The study builds on the method 
of analysis developed in a previous paper for one source, GL 2591, by van dcr Tak et al. (1999). 
The column density is inferred from continuum observations, and the density from molecular lines, 
while the temperature is calculated self-consistently from the luminosity, dust properties, and n(r). 
We use the derived density and temperature profiles to characterize the early evolution of massive 
stars and their surroundings. 

2. Choice of Targets 

The fourteen sources studied in this paper are introduced in Table 1, which lists the source 
names, positions, distances, luminosities, radio continuum flux densities and associated IRAS Point 
Source Catalog (PSC) entries. All sources have been well studied before at radio and infrared wave- 
lengths. They were selected to be luminous (> 10^ Lq) and visible from the Northern hemisphere. 
The distances to some of the sources are quite uncertain; in particular, for the sources in the 
Cygnus X region, we use a fiducial distance of 1 kpc (van der Tak et al. 1999), for ease of scaling 
once better distances are available. 

The main sample consists of nine deeply embedded massive young stars, which were addi- 
tionally selected for mid- infrared brightness (> 100 Jy at A12 fim), allowing comparison to exist- 
ing absorption line studies. A combined analysis of emission and absorption lines gives powerful 
constraints on the physical structure and geometry of the system, and allows a nearly complete 
inventory of the chemical composition of both the solid and the gas phases. These nine sources 
have been extensively studied in the infrared, both from the ground (Willner et al. 1982; Mitchell 
et al. 1990) and with the Infrared Space Observatory (van Dishoeck et al. 1998). As a comparison 
sample, five luminous embedded young stars which are weak in the mid-infrared are studied, whose 
surroundings contain a "hot core" and/or an ultracompact H II region. 

The full sample contains sources with luminosities between 1x10^ and 2 x 10^ Lq and distances 
between 0.9 and 4 kpc. The radio flux densities are at the ~mJy level for most sources, which is 
orders of magnitude lower than expected from an H II region photoionized by a star with the same 
luminosity. The radio emission from the bright mid-infrared sources has a spectral index 7 in the 
range 0.5 — 1.0 and may arise in a spherical wind for which 7 = 0.7. The principal difference between 
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our main sample and that of other studies of massive young stars (e.g., Tofani et al. 1995; Hunter 
1998; Sridharan et al. 1999) is the selection on mid-infrared brightness, which could introduce an 
orientation bias, in the sense that low-density cavities such as produced by molecular outflows 
may be preferentially directed towards us. Alternatively, the high mid-infrared brightness may be 
an evolutionary effect. Depending on the relative importance of different effects, the mid-infrared 
bright sources in the main sample could be younger or older than the sources in the comparison 
sample. We will discuss this point in more detail in Section 6.5. 

The sources were selected on isolated location within ~ 30" at infrared wavelengths, so that 
the heating is dominated by the central star. However, while GL 2591 and GL 490 are examples 
of objects forming in isolation, NGC 6334, S 140, NGC 7538 and W 3 are regions where low-mass 
stars are also born, as revealed by near-infrared images (Straw & Hyland 1989; Evans et al. 1989; 
Hodapp 1994; Megeath et al. 1996; Bloomer et al. 1998). In addition, S 140 contains three massive 
objects separated by 10 — 15", which contribute about equally to the far-infrared luminosity (Evans 
et al. 1989). For the other sources, this information is not available. In the case of W 3 IRS5, 
radio continuum observations (Claussen et al. 1994) suggest the presence of seven early B-type 
stars within 4", a very large number in view of the average mass function and stellar density in 
young clusters (Carpenter et al. 1993; Luhman et al. 1998; Alvcs et al. 1998). As long as the power 
source is much smaller than our beam, its multiplicity does not matter for our analysis. 

In several cases, high-resolution continuum observations reveal the presence of an (ultra- 
) compact H II region next to the source of interest. Since luminosity and/or mass are dominated 
by the infrared source, our models of centrally heated sources are still reasonable approximations. 

3. Observations 

3.1. Single Dish Submillimeter Spectroscopy 

Spectroscopy of selected molecular lines in the 230, 345 and 490 GHz bands was performed 
with the 15-m James Clerk Maxwell Telescope (JCMT)^ on Mauna Kea, Hawaii during various runs 
in 1995-1997. In total, 15-20 lines were observed per source. The antenna has an approximately 
Gaussian main beam of FWHM 18" at 230 GHz, 14" at 345 GHz, and 11" at 490 GHz. Detailed 
technical information about the JCMT and its receivers and spectrometer can be found on-line at 
<http : //www . j ach . hawaii . edu/ JCMT/home . html>. Receivers A2, B3i and C2 were used as front 
ends at 230, 345 and 490 GHz, respectively. The Digital Autocorrelation Spectrometer served as 
backend, with continuous calibration and natural weighting employed. Values for the main beam 
efficiency, rj^b, determined by the JCMT staff from observations of Mars and Jupiter, were 0.69, 



^The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre, on behalf of the Particle Physics 
and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Scientific Research and 
the National Research Council of Canada. 



-6- 



0.58 and 0.53 at 230, 345 and 490 GHz for the 1995 data, and 0.64, 0.60 and 0.53 for 1996 and 
1997. Absolute calibration should be correct to 30%, except for data in the 230 GHz band from May 
1996, which have an uncertainty of ~ 50% due to technical problems with receiver A2. Pointing 
was checked every 2 hours during the observing and was always found to be within 5". Integration 
times are 30-40 minutes per frequency setting, resulting in rms noise levels in Tmb per 625 kHz 
channel ranging from ^ 30 mK at 230 GHz to 100 mK at 490 GHz. 

To subtract the atmospheric and instrumental background, a reference position 180" East was 
observed using the chopping secondary mirror. For the C^^O J = 2 ^ 1 line, we also position- 
switched using an 1800" offset, which increased the line fluxes typically by 15%. The contribution 
by extended emission in other lines should be less since they all need higher densities for the 
excitation. Where two measurements of C^^O J = 2 — > 1 are available, the position-switched data 
will be used. 

Observations of the C'^^S J = 2 — > 1 and J = 3 ^ 2 and CS J = 5 ^ 4 lines were carried out 
with the 30-m telescope of the Institut de Radio Astronomic Millimetrique (IRAM) on Pico Veleta, 
Spain, on January 28-30, 1999. Receivers BlOO, C150 and B230 were used at 96, 145 and 245 GHz, 
respectively, tuned single sideband. The autocorrelator was used as the backend, which covers a 
bandwidth of 110 - 170 km s""^ at a resolution of 0.2 km s'^ The FWHM beam sizes are 24", 16" 
and 10.4"; the forward efficiencies were 0.90, 0.82 and 0.84 and the beam efficiencies 0.75, 0.55 and 
0.48. The data were calibrated onto the Tmb scale by multiplying T| by the ratio of the forward 
and beam efficiencies. As a calibration check, observations of GL 2591 and NGC 7538 IRSl were 
compared to results from Plume et al. (1997) and found to agree to 10%. 

The data on GL 2591 have been presented in van der Tak et al. (1999). The observations on 
GL 490 are given in Schreyer et al. (2000, in prep.). Most of the data on W 3 (H2O) and W 3 IRS5 
are taken from Helmich & van Dishoeck (1997). In addition, we used data from the surveys by 
Anglada et al. (1996), Bronfman et al. (1996), Plume et al. (1992, 1997), and from the observations 
of Kastner et al. (1994) for GL 2136, Zhou et al. (1994) for S140 IRSl, Cesaroni et al. (1997) for 
IRAS 20126, Hauschildt et al. (1993) and Mangum & Wootten (1993) for DR 21 (OH), Dartois 
(1998, Ph.D. thesis) for GL 7009S and Olmi k Cesaroni (1999, A&A, in press) for W 28 A2. Care 
was taken not to include data at positions > 5" away from those in Table 1. For the sources 
GL 7009S, IRAS 20126+4104 and W 28 A2, no H2CO data are available. 

Line parameters were measured by fitting a Gaussian profile, where the free parameters were 
the total line flux J TuBdV, the FWHM line width and the center velocity. The measured line 
fluxes are collected in Table 2; Table 3 gives for each source the central LSR velocity and the 
FWHM line width. The values are the averages of the C^^O and C^^S lines (4 lines in most cases), 
which are assumed to be optically thin. The error bars reflect the spread among the lines, which 
for all sources exceeds the error in the individual measurements for both quantities. None of the 
CS and H2CO line profiles show obvious evidence of self-absorption. 
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3.2. Single-dish Mapping 

Maps of the CS J = 5 ^ 4 and/or 7^6 and/or C^^O J = 2 ^ 1 lines were obtained in 

1998 July and December and 1999 July with the 10.4-m antenna of the Caltech Submillimeter 
Observatory (CSO)^. The backends were the facility 500 and 50 MHz bandwidth Acousto-Optical 
Spectrometers (AOS). The pointing is accurate to within 4". Beam sizes and main beam efficiencies 
were 32" and 0.66 at 230 GHz and 21" and 0.61 at 345 GHz. 

In July 1998, submillimeter continuum maps at A350 fim of W 33A, GL 2136, S 140, GL 490, 
GL 2591 and W 3 IRS5 were made using the Submillimeter High Angular Resolution Camera 
(SHARC) (Hunter et al. 1996) on the CSO. The CSO has a beam of size 10" at this wavelength. 
The weather was good, with zenith optical depths at A350 /nm of 0.92 — 1.66. From observations 
of Uranus, the gain was measured to be (140 ± 30) Jy/V. 

3.3. Interferometer Observations 

Maps at 86-230 GHz of GL 2136, NGC 7538 IRSl, NGC 7538 IRS9 and W 33A were obtained 
with the six-element interferometer of the Owens Valley Radio Observatory (OVRO)^. The OVRO 
interferometer consists of six 10.4 m antennas on North-South and East- West baselines. A detailed 
technical description of the instrument is given in Padin et al. (1991). Data were collected in 1997- 

1999 in several array configurations at frequencies of 86, 106 and 115 GHz. The two sources in 
NGC 7538 were observed in the compact and extended array configurations, while for the Southern 
sources GL 2136 and W 33A, a hybrid configuration with long North-South and short East- West 
baselines was also used to improve the beam shape. Baseline lengths range from the shadowing 
limit out to 80 kA at 86 GHz and out to 150 kA at 230 GHz, corresponding to spatial frequencies 
of -2500 to ~ 10^ rad-^ 

The antenna gains and phases were monitored with short integrations of nearby quasars: 

BL Lac for NGC 7538 IRSl and IRS9, and NRAO 530 for W 33A and GL 2136. Passband 
corrections were derived using a local signal generator and by fitting first order polynomials to data 
on 3C273. Flux calibration is based on snapshots of Uranus and Neptune, and is estimated to be 
accurate to 30 % at 86 - 115 GHz and to 50 % at 230 GHz. 

Simultaneous continuum observations in the 230 GHz window were also performed, but only 
in the cases of W 33 A and NGC 7538 IRSl did the weather allow useable data to be taken. In 
addition, observations of molecular lines at 86 — 230 GHz were carried out with the OVRO digital 



^The Caltech Submillimeter Observatory is operated by the California Institute of Technology under funding from 
the U.S. National Science Foundation (AST96-15025) 

^The Owens Valley Millimeter Array is operated by the California Institute of Technology under funding from the 
U.S. National Science Foundation (AST96-13717). 
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correlator, the results of which will be presented elsewhere. 

4. Results 

4.1. Molecular Emission Line Profiles 

The CS line profiles are presented in Figure 1. In addition to a strong, single-peaked line core, 
which is also detected in isotopic lines (C^^S, ^^CO) and which has an approximately Gaussian 
shape, high-velocity wings are clearly detected. These wings are more prominent in J = 5 — > 4 

than in highcr-J lines, are also detected in the strongest H2CO lines, and must arise on small 
scales, since they are stronger in the IRAM 30m beam than in the JCMT beam. The wings are not 
detected in the W 3 sources, for which we do not have J = 5 ^ 4 data. Since the sources studied 
in this paper are known to drive CO outflows, it is natural to associate the CS and H2CO wings 
with dense gas in these outflows. 

In all bright mid-infrared sources studied in this paper except W 33A, blueshifted infrared 
absorption of CO is detected at similar velocities as in CO rotational line emission (Mitchell et al. 
1991), and in the cases of W 3 IRS5 and GL 2591 up to much higher velocities, implying that the 
outflow is directed along the line of sight. In realistic outflow models with a finite opening angle, 
a mixture of blue- and redshifted emission is expected on both sides of the driving source, which 
should be visible in submillimeter emission. However, the wings seen on the CS and H2CO lines are 
much stronger at blueshifted than at redshifted velocities, and for some sources no redshifted wing 
emission is detected at all. Since the line asymmetry is stronger in the IRAM 30m spectrum of CS 
J = 5 ^ 4 than in the JCMT profile of the same line, the asymmetry must arise on scales of ^ 10". 
Since the redshifted outflow lobe lies in the background, we suggest that the lack of redshifted 
CO and CS wing emission is due to obscuration. Since the outflow lies at a velocity offset, the 
absorption must be continuous absorption by dust. The H2 column density required to do this is 
~ 10^^ cm~^ on a J; 10" scale, corresponding to a visual extinction of ~ 10^ magnitudes. 

4.2. Average Physical Conditions 

Estimates of the mean temperature and density of the gas can be obtained by comparing 
observed line ratios of specific molecules with non-LTE models which include radiative trapping. 
Examples of this approach can be found in Jansen et al. (1994) and Helmich et al. (1994); it is 

similar to the "Large Velocity Gradient" method used by Zhou et al. (1994) and Plume et al. 
(1997). The observed line ratios of C^^S are well-suited to constrain the H2 density. We have 
calculated synthetic line ratios for a range of temperatures and densities, using the rate coefficients 
by S. Green (priv. comm., cf. Turner et al. 1992). Comparison with the data indicates densities of 
10^ to > 10^ cm~^, but different line ratios generally give inconsistent answers for the same source. 
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indicating the presence of density inhomogeneities, probably in the form of a gradient since hnes 
observed with several telescopes are systematically brighter when observed with a smaller beam. 

The measured ratios of the H2CO J = 624 ^ 423/J = 5o5 4o4 and the J = 642/41 ~^ 
^41/40/ J = ^24 lines are particularly sensitive to temperature. Unlike C^^S, these H2CO hne 

ratios were measured in the same beam, or even in the same spectrum, improving the calibration 
between pairs of lines, although their filling factors might still be different. Model calculations 
using collisional rate coefficients by Green (1991) give temperatures of 60 K for S 140 IRSl to 
> 200 K for W 3 (H2O) and NGC 6334 IRSl. Again, the two ratios usually do not yield the 
same temperature for the same source, indicating the need for models with a varying temperature. 
In general, the modeling indicates somewhat higher temperatures and densities for the weak mid- 
infrared sources. In Section 5, we will see that this effect is due to a steeper density gradient in 
these sources. 

The CO column densities derived from C^^O, using the temperatures and densities found above, 
are listed in column 3 of Table 3. A plot of these values against the column densities derived by 
Mitchell et al. (1990) from infrared absorption observations is shown in Figure 2. Abundance ratios 
of i^C/i^c = 60 and ^^o/i^O = 2500 (Wilson & Rood 1994) are assumed. The two measurements 
agree to a factor of 2 for all sources, and often much better. This result implies that the circumstellar 
envelopes of our objects have to first order a spherical shape, with the infrared source in the center 
(see also Section 6.1). The average absorption column density is 66% of the emission value, which is 
higher than half as expected in a uniform model. Beam dilution in the emission data could explain 
this difference, in which case the sources are centrally condensed. 



4.3. Submillimeter Continuum and Line Maps 

The CS and C^^O maps are presented in Figure 3; Figure 4 shows the SHARC A350 /xm 
maps and maps of NGC 7538 at AA450, 850 /xm, obtained with the Submillimeter Common-User 
Bolometer Array (SCUBA) on the JCMT, provided by G. Sandell (1998, private communication). 

The maps appear compact but slightly extended. In a few cases, the map peak is offset from 
the center position, but this offset falls within the 5" pointing uncertainty. In the case of CS, this 
suggests that the infrared sources coincide with local density maxima in the surrounding molecular 
cloud. The C^^O maps are sometimes peaked at the same position as CS, which implies a maximum 
of the column density, but more often, these maps are not strongly peaked and have a much lower 
dynamic range than the CS maps. The map of NGC 7538 shows a chaotic structure, rather than 
clear peaks. 

Column 4 of Table 3 lists the diameter of the CS J = 5 — > 4 emission, measured from the 
maps as the point where the brightness has dropped to 50% of the maximum. These numbers will 
be used in the next section to constrain the radii of the models. The values span a fairly narrow 
range: 36 — 64" and do not show an obvious correlation with distance, suggesting that the cores 
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do not have a single intrinsic linear size. Instead, it appears that the molecular gas has a scale-free 

density structure, such as can be described by a power law. This conclusion is supported by the 
fact that the CS diameters are somewhat (10 — 100%) larger than the beam size, which is also what 
power law models predict. 

The maps do not show a cutoff in the emission, such as would be produced if the star-forming 
cores had a distinct edge. Rather, the brightness keeps going down until the detection limit. Maps 
of a wider field and with a higher dynamic range than presented here may reveal if there is such an 
edge or if at large radii, the density gradient flattens out and the core merges into a surrounding 
molecular cloud. Some objects, such as W 3 IRS5, DR 21 (OH), NGC 6334 IRSl and the NGC 7538 
sources, are clearly part of a more extended giant molecular cloud. However, the focus of this paper 
is on the ~ 30" region around the massive young stars. 

4.4. Interferometric Continuum Maps 

Figure 5 presents maps created from the OVRO data by gridding and Fourier transforming 
the visibility data with uniform weighting. Deconvolution with the CLEAN algorithm and self- 
calibration of the uv phases based on the brightest CLEAN components helped to improve the 
image quality. Table 4 lists the positions and flux densities of the detected sources. 

The maps typically show a single compact source at or near phase center, within the ~ 5" 
positional error of the infrared positions from Table 1. The extended 107 GHz emission North of 
NGC 7538 IRSl is the H II region surrounding NGC 7538 IRS2, which has also been detected at 
lower frequencies with the VLA (e.g., Henkel et al. 1984). In the case of W 33A, two sources are 
detected, neither of which coincides with the IR position of Dyck & Simon (1977). The brightest 
source, MMl, coincides however with a VLA detection (Rengarajan & Ho 1996) and with the 
infrared position by Capps et al. (1978), who also report a second 2.2 /xm source 3" South of W 33A, 
which may be related to MM2. It is unknown if these two sources are physically associated; if they 
are, they may represent a young massive binary star at a separation of 4'.'7 or 19, 000 AU. This 
separation is larger than the ~ 10^ AU found by Wyrowski et al. (1999) for W 3 (H2O) and by 
Mundy et al. (1999) for a sample of low- mass young stellar objects, but comparable to that found 
by Woody et al. (1989) and Padin et al. (1989) for DR 21 (OH). 

The flux densities found here for NGC 7538 IRSl are in agreement with those by Woody et al. 
(1989) and Akabane et al. (1992). The spectral index, measured from 107 to 230 GHz, is (0.9 ±0.5). 
The spectral indices of the sources in W 33A are (1.8 ± 0.6) (MMl) and (1.6 ± 0.6) (MM2). The 
values of the spectral indices for W 33A rule out optically thin emission from either ionized gas or 
from dust. They are, however, consistent with black body emission, probably from a compact dust 
source. In Paper 1, OVRO observations of GL 2591 were presented which gave similar results; we 
refer the reader to that paper for a detailed discussion. 
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5. Models 

Motivated by the results of the excitation analysis (§ 4.2) and of the CS maps (§ 4.3), we will 
proceed by developing spherically symmetric power law models. The modeling procedure follows 
van der Tak et al. (1999) to which we refer for details. 

5.1. Dust Continuum Models: Mass and Temperature Structure 

The dust emission from the sources was modeled using the one-dimensional diffusion code by 

Egan et al. (1988). A density structure of the form n = no(r/ro)~" was used, with a in the range 
0.5 — 2.0. The fiducial radius ro was set to be half the diameter of the CS emission {D(CS) in 
Table 3), converted to AU in Table 5. The density no at this radius was derived for each a by 
matching the submillimeter dust emission; in the next section, we will constrain the value of a by 
modeling molecular lines. To avoid edge effects in modeling the full range of emission, we used 
outer radii in our models that are twice tq. The inner radius was arbitrarily set to 1/300 of the 
outer radius, small enough not to influence the calculated brightness or temperature profile, as 
verified by test calculations for W 33A and GL 2591. Dust opacities appropriate for star-forming 
regions were taken from Osscnkopf &: Hcnning (1994) (their Model 5) and converted to absorption 
cross-sections using a grain mass density of 2.5 g cm^'^ and a radius of 0.1 /xm. This radius is 
close to the median of more realistic size distributions, so that the calculated dust temperatures 
represent the bulk of the dust (cf. Churchwell et al. 1990). Dust properties are known to vary 
from one region to the next (Lis et al. 1998; Visscr ct al. 1998; Hogcrhcijdc & Sandell 2000) and 
are likely to change inside our sources as well, due to the changing temperature (MenncUa ct al. 
1998; see § 6.2). Maps at several far-infrared wavelengths at ^10" resolution would allow us to 
disentangle temperature and dust opacity variations, but awaiting such data, we assume a constant 
grain opacity. With the luminosity, the distance and the size of the source fixed, M(< ro) was used 
as the only free parameter to match the observed submillimeter continuum emission. 

Figure 6 shows the synthetic continuum spectra compared with observations. Most submil- 
limeter data are from this work and from Sandell (1994); the sources of the additional data are 
listed in the caption. The models have been constructed to fit the submillimeter (A300 — 1300;um) 
data, and are seen to reproduce the observed emission at Ai^50 fim for every source, i.e. up to 
optical depths of a few. The shorter-wavelength emission is in general not matched, although care 
was taken to include only photometry in small (^ 10") beams to avoid a contribution from reflection 
nebulosities. The high brightness of GL 2591 at A20/Ltm was attributed by van der Tak et al. (1999) 
to the evaporation of ice mantles close to the star, which decreases the A20/im optical depth by 
30% (Osscnkopf & Henning 1994). Similar effects play a role for the sources presented here, as 
illustrated for W 33A in the figure. The emission at A$ lO^um is not expected to be well reproduced 
in these models, since the high optical depth makes it very sensitive to deviations from the assumed 
spherical shape. The fact that all sources require less extinction at short wavelengths is discussed 
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quantitatively in § 6.1. This point is illustrated by the fact that the column densities derived in this 
paper are much higher than those found by Faison et al. (1998) by fitting the near- to far- infrared 
spectra of a sample of ultracompact H II regions, which have similar submillimeter flux densities 
and lie at similar distances. 

The calculated temperatures follow an r~^''^ profile at distances greater than 2 — 3 x 10^ AU 
from the star, with the absolute temperature scale set by the luminosity. Inside this radius, the 
envelope is optically thick to the photons carrying most of the energy, and the temperature gradient 
is steeper than r~°-^, with the absolute scale set by the extinction, which acts as a blanket. Hence, 
for a given luminosity, a lower extinction leads to a higher temperature at a certain radius. This 
is the case for our models with shallow density gradients. For a given column density, a shallower 
density gradient implies a lower extinction (in a pencil beam), because a larger fraction of the beam 
is filled with warm dust. We therefore expect, like Chini et al. (1986) and Churchwell et al. (1990), 
that bright mid-infrared sources have shallow density gradients. However, it is seen from Figure 6 
that even the models with a = 0.5 fail to fit the near-infrared part of the spectrum. This is not 
due to our selection on brightness in the mid-infrared since the data of the comparison samples are 
not matched either. The reason for the discrepancy must be sought in (small) geometrical effects. 



To determine the slope of the density gradient a, the C O, CS, C S and H2CO line radiation 
from the 14 sources was modeled with the Monte Carlo code written by Hogerheijdc & van der Tak 
(2000). The models consist of 30 spherical shells spaced logarithmically within the same radii as the 
dust models. The same density structures were used as in the dust models, and the temperature 
structure was taken from the dust continuum models for the same value of a. The intrinsic line 
profile was assumed to be a Gaussian with the measured FWHM (Table 3). Initially, molecular 
abundances were assumed to be constant throughout the model. 

The gas column density (or equivalently M(< ro) or the value of no) for each a was also taken 
from the dust models, assuming a gas to dust mass ratio of 100. This assumption was tested for 
GL 2591 against C^^O observations by van der Tak et al. (1999) because this source has negligible 
depletion of CO (see also § 6.2). With the column density fixed, the density profile can be obtained 
by modeling the emission lines of the high-critical density molecules CS and H2CO. After solving 
for the molecular excitation, velocity-integrated intensities have been calculated in the appropriate 
beam for each observation. In many cases, a line was observed in several beams, which in our models 
acts as a simple substitute for simulating map data. Comparison to the observations proceeds by 
minimizing the quantity 



where F^"^''^ and ^(modei) 

are the observed and synthetic line fluxes and the sum is taken over all 



5.2. Models of the Line Emission: Density Structure 
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N observed lines. A grid of models was run for each molecule by varying the density law exponent 
a and the molecular abundance. 

The success of any fit procedure depends critically on a good understanding of the error budget 
((t). Most line fluxes were assigned a 30% error, but the following lines have a larger uncertainty 
of 50%: (i) lines in the 460-490 GHz band, because of the more uncertain calibration; (ii) lines 
measured in beams of ^30" diameter, which may suffer from confusion; and (iii) the J = 1 — > 
lines of CS and C^^S, which may contain a contribution from the surrounding cloud. The following 
observations were found to be more than a factor of 2 off any of our models, and therefore discarded 
in the calculation: CS J = 5 ^ 4 and 7 ^ 6 in GL 490 and CS 7 ^ 6 in W 33A from Plume 
et al. (1997), CS 5 ^ 4 and C^^S 3 ^ 2 in S140 IRSl from Zhou ct al. (1994), C-'^^S 2 ^ 1 in IRAS 
20126+4104 from Cesaroni et al. (1997) and C^^S J = 7 ^ 6 in W 28 A2 from Plume et al. (1997). 

The results of the emission line models are presented in Figure 7. The fit quality parameter 
is plotted as a function of the density law exponent a and the CS abundance with respect to 
H2. The elongation of the contours indicates that the quality of the fit depends to first order on 
the CS abundance and to second order on the density profile. The parameters of the best fitting 
models are summarized in Table 5. 

It is seen from Table 5 that for the main sample, a = 1.0 — 1.5, while the other sources have 
a = 1.75 — 2.0. This result is consistent with the result from the dust models (§ 5.1), namely 
that for a given source, the model with lowest a has the highest mid-infrared flux. Thus, the 
mid-infrared brightness of the main sample is not a pure orientation effect. However, since the 
spherical dust models fail to reproduce the ncar-infrarcd emission of all sources, deviations from 
spherical symmetry are important. This conclusion is supported by synthetic CS line profiles from 
our power law models, which are self-absorbed for all but the least massive sources. The effect of 
deviations from spherical symmetry on line profiles and near-infrared emission was discussed for 
GL 2591 in van der Tak et al. (1999). Decreasing the central column density by a factor of a few 
can increase the mid-infrared emission and decrease the self-absorption substantially. 

The best-fit CS abundance is 3 x 10~^ on average, with source-to-source variations of a factor 

of 3 for most sources. The hot cores and the ultracompact H II region have the highest CS 
abundances: (1 — 2) x 10~^. The low CS abundance in GL 7009S, 0.4 x 10^^, is probably due to 
freezing out of the molecules onto the dust grains, as STiggcstcd by the large column densities of 
several ice species towards this source (d'Hendecourt ct al. 1996). 

Figure 7 also lists the number of CS and C'^'^S lines observed. The sources for which the most 
lines have been observed also have the tightest constraints on a. For less well-observed sources, 
our results may be influenced by sampling effects. In particular, if only a small range of molecular 
energy levels is available, the value of a may be biased, or a gradient may be hard to detect. We 
have checked for such effects by calculating J, the average upper J-level of the data set. Most 
sources have J = 4.5 — 5.5, which implies that the observable range of critical densities has been 
evenly sampled. However, for GL 7009S, and IRAS 20126, for which fewer high-excitation line data 
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exist, J = 3.2 — 3.6, and the model results axe therefore less robust. 

The value a = 2 found for the hot cores is higher than for the other sources, and we have 
investigated if this result could be an artifact of our assumption of a constant molecular abundance. 
Hot cores have high abundances of many molecules including saturated organics, which are thought 
to be due to freshly evaporated grain mantles. In this picture, the radius where the temperature 
reaches 90 K is an important threshold because water ice, the main ingredient of interstellar grain 
mantles, evaporates. It is possible that the abundance of CS is enhanced above this temperature as 
well. To test this idea, a model has been run with the luminosity and the M(< ro) of W 3 (H2O), 
but with a = 1.5. The CS abundance is 5 x 10~^ in the outer parts, but enhanced where the 
dust temperature exceeds 90 K. A match to the data of equal quality as with the best constant- 
abundance model (a = 2.0, CS/II2 = 1 x 10~^: Table 5) was obtained by enhancing the CS 
abundance in the inner region by a factor of 10. This result indicates that the hot cores may have 
a density structure similar to that of the main sample, if plausible variations of CS/H2 with radius 
are considered. SCUBA maps of other "hot core"-type sources by Hatchell et al. (2000, submitted 
to A&A) indeed suggest values of a f« 1.5. 



5.3. Comparison of Dust and Gas Tracers 

The density structures derived from the CS excitation will now be tested by using them to 
model the radial profiles of the submillimeter dust emission. This emission is optically thin and 
hence measures the column density distribution, given the calculated temperature structure. The 
points in Figure 8 are averages of slices along the North-South and East- West directions through 
the images shown in Figure 4. To reduce the noise, the data have been folded about the image 
maximum. In the case of W 3 IRS5, the local maximum in the SHARC map that corresponds to 
the infrared position was used to center the slices, and the Western direction was not included in 
the scans because of confusing extended emission clearly visible in Figure 4. Superposed are slices 
through model images for various values of a as dotted lines, while the solid line corresponds to 
the value of a derived in § 5.2 from the CS excitation. The model profiles (in a 10" beam) were 
calculated with a code kindly supplied by L. Mundy. The models are exactly those used to produce 
Fig. 7, except that for W 3 IRS5, GL 490 and NGC 7538 IRS9, the outer radii were increased to 
60" to avoid edge effects. In these cases, the radial profile of the dust emission was quite extended. 

Figure 8 demonstrates that the dust and gas tracers agree very well on the best value of a, 
implying that the volume density and column density distributions are consistent. This result 
indicates that the dust and gas are well-mixed and that the structure of the envelopes is fairly 
homogeneous and not very clumpy. The only exception to this rule is S 140 IRSl, for which the 
dust gives a = 1.0 while the gas gives a = 1.5. Although an inhomogeneous (clumpy) structure 
would cause a discrepancy in this direction, it is hard to see why only one of our sources would 
have such a structure. Clumps in S 140 have been proposed by e.g. Spaans &; van Dishoeck (1997), 
but in a much larger (~arcmin-sized) region outside the dense star-forming core. More likely to be 
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important here is an elevated temperature caused by the nearby sources IRS2 and IRS3 at 10 — 15" 

offsets, which have luminosities similar to that of IRSl (Evans et al. 1989), and by external heating 
by the nearby ionization front (Lester et al. 1986). So S 140 is unique in our sample in that it is 
the only source where our assumption of one central heating source breaks down. 

Although the data on GL 2591 arc well matched by the power law on average, the slope of 
the data is shallower inside a radius of 10 — 20", and drop more steeply outside this radius. This 
suggests a variation in the temperature or the density gradient with radius. Since the effect is 
more pronounced towards shorter wavelengths, a variation of the temperature gradient with radius 
seems more plausible. An inner region of roughly constant temperature, such as caused by a central 
cavity with little or no extinction, may reproduce the data, but such models fall outside the scope 
of this paper. 

6. Discussion 

We have modeled single-dish continuum and line data with power law models n = nQ{r/rQ)~", 
and obtained good solutions (x^ = 1 — 3) with a ~ 1.0 — 1.5 for the bright mid-infrared sources and 
a = 1.5—2.0 for the other sources, although the highest value a = 2.0 found for two "hot cores" may 
be an overestimate due to an enhanced CS abundance close to the center. This section investigates 
the validity of the assumptions that went into the models and discusses possible implications of 
their results. 



6.1. Envelope Structure 

The models developed in Section 5 assume that the clouds arc spherical and homogcncoTis. 
Mitcheh et al. (1990) observed the ^^CO v = 1 ^ band at 4.7 yUm in absorption towards the 
bright mid-infrared sources. The excitation of the quiescent (non-outflowing) gas rules out single- 
temperature models, but the data are well fitted by a model with two discrete temperatures. The 
column density ratio of these "cold" and "hot" components varies from 1 to !=a 5. The physical 
origin of such a two-temperature structure is not clear, however. We have calculated the column 
densities in each rotational state of ^'^CO up to J = 25 from our power-law models, and the results 
are compared to the data of Mitchell et al. in Figure 9. The low-resolution data for GL 7009S from 
Dartois et al. (1998b) may contain a contribution from outflowing gas, and are hence regarded as 
upper limits. The models have been scaled to the observed total column density, which is always 
within a factor of 2 from the column density measured by our submillimeter emission data (§ 4.2). 
The good match (maximum deviation of a factor of 3) to the infrared data, which sample a pencil 
beam toward the infrared source, indicates that a spherically symmetric model is a reasonable first- 
order description of these sources. Together with the evidence from the CS emission line profiles 
and the near-infrared continuum emission, we estimate the deviations from spherical geometry as 
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a decrease in optical depth by a factor of 3 over an < 10" area. 

The assumption of a homogeneous density structure is tested in Figure 10, where M(< ro) is 
compared to the mass derived from the virial theorem inside the same radius. The latter method 

assumes that the cloud's gravitation just balances its pressure as measured by the line width AV 
(Table 3). If the gas is clumpy, the virial mass is smaller than the power law mass by a fraction f^, 
with the volume filling factor < f\, < 1. The ratio Mv/M{< ro) has a mean value of 2.77 and a 
standard deviation of 1.66. The only source for which M(< ro)> My is GL 7009S, for which the 
model results are uncertain. 

While there is no evidence for clumping within the envelopes, other authors have found that 
regions of massive star formation are usually clumped on a much larger scale. It may not be a 
coincidence that the sizes derived for the clumps in S 140 by Plume et al. (1994) and in M 17 SW 
by Wang et al. (1993), 0.2 pc, are similar to the sizes of the envelopes found here. It is however 
outside the scope of this paper to investigate a possible evolutionary connection between molecular 
cloud clumps and the envelopes of massive young stars. 

6.2. Relation of CO Abundance with Temperature 

Our assumed values of the gas to dust mass ratio and the dust opacity at submillimeter 
wavelengths were tested by modeling C^^O J = 3 — > 2 and 2^1 observations. The best-fit values 
for the abundance of CO are listed in column 6 of Table 5. The derived abundances are 10 — 60% 
of the value of 2.7 x 10~^ measured by Lacy et al. (1994) towards the warm cloud NGC 2024. Our 
CO/H2 value for GL 2591 is in good agreement with the value of 3.1 x 10~^ measured by C. Kulesa 
(1999, priv. comm.), while in the case of GL 490, our value is a factor of 7 lower. The CO abundance 
measured in emission is expected to be lower than in absorption, because the emission data are 
more sensitive to extended cold gas. When measured in infrared absorption, which samples warm 
material close to the star, the column density ratio of solid to gaseous CO never exceeds unity 
towards massive young stars (Mitchell et al. 1990; van Dishoeck 1998). 

Although source-to-source variations in grain opacity cannot be ruled out, a more likely ex- 
planation of the spread in CO abundances is that different amounts of CO are frozen out on grain 
surfaces. This occurs at temperatures ^20 K in the case of pure CO (Sandford & Allamandola 
1993), and up to 45 K for CO-H2O mixtures. In both correlation of the abundance 

with temperature is expected. Figure 11 shows a plot of the CO abundance derived from C^''0 
observations versus the mass-weighted temperature, defined as T = J T {r)n{r)r^ dr / J n{r)r^dr (see 
Table 5). The two quantities are seen to be correlated, which we interpret as the effect of depletion 
of CO at low temperatures. 

These results strongly suggest that depletion and thermal desorption are the processes con- 
trolling the abundance of CO in the gas phase, after chemical reactions lock up almost all of the 
available carbon. This situation is in contrast with that of CS, for which the abundance shows no 
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clear trend with T. This difference presumably reflects the much greater chemical inertness of CO 
compared to CS, and the lower evaporation temperature of CO compared with CS. 

In the case of GL 7009S, the ^^CO column density observed in absorption by Dartois et al. 
(1998b) is similar to that in our model. This implies that along this line of sight, CO is much less 
depleted than CS, for which a depletion by a factor ~ 10 was found compared to the other sources 
in our sample (Section 5.2). Again, this can be understood from the difference in evaporation 
temperature for these two molecules. The large column of CO ice towards, e.g., W 3 IRS5, is 
not predicted by our model. We suggest that this star-forming core is surrounded by an extended 
dense cloud which also produces the self-reversed CO emission proflles. Sources like GL 2591 and 
GL 2136 do not seem to have such a "skin" , or they have at most a translucent one. 

Mennella et al. (1998) found that the submillimeter opacity of dust grains increases by 10—50% 
when the temperature rises from 24 to 300 K. This effect would decrease the masses of our sources 
and increase the inferred abundances of CO, and qualitatively produce the same trend as seen in 
Figure 11. However, over the applicable temperature range, 20 — 50 K, an effect of only 5 — 25% 
should occur, which is much less than the factor of 3 increase in CO abundance that we flnd. 
Although the grain opacity is likely to vary within our sources, they do not influence the conclusion 
that the CO abundance is controlled by freezing and sublimation. 

6.3. The Inner < 1000 AU: Evidence for Compact Dense Material 

In this section, the envelope models derived in § 5 will be compared to the OVRO continuum 
data. Since the images presented in Figure 5 show only compact, circularly symmetric structure, a 
simple one-dimensional analysis in the Fourier plane is sufficient. 

The points in Figure 12 are the OVRO data, averaged in annuli around the source in the uv 
plane. In the case of W 33A, the source MM2 was subtracted from the data before averaging. The 
error bars represent the la spread between the data points in each bin, and do not include the 
overall calibration error. Superposed on the data points are model curves, derived by calculating 
the dust emission from the best-fit models (Table 5) and Fourier transforming the result. 

Only at the lowest spatial frequencies, the data match the model curves within the calibration 
error. Most data points lie well above the model curves, and the observed amplitudes do not drop 
with increasing spatial frequency, as the models do. These differences suggest that the emission 
detected with OVRO is not related to the envelopes seen in the single-dish data, but due to compact 
structure of size $2". The same conclusion was reached for GL 2591 by van der Tak et al. (1999). 

Towards the highest observed spatial frequencies, 0.5 — 1 x 10"'' rad~^, corresponding to ~ 
3000 AU at 1 kpc, the compact emission is 10 — 100 times stronger than that of the envelope model. 
This factor is a lower limit to the column density contrast between the two components, since the 
compact emission is probably optically thick, as indicated by the spectral indices. Combined with 
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the envelope column densities in Table 3, this result suggests that the compact emission seen with 
OVRO can explain the asymmetry observed in the CS line profiles (§ 4.1). 

Further constraints on the nature of the compact emission can be obtained from the flux 
densities (Tabic 4), which imply brightness temperatures of 50 — 80 K for NGC 7538 IRSl and 
« 1 K for the other sources. These values imply that the emission is beam diluted by a factor 
of ^ 100, since the physical temperature of the compact structure must be at least that of the 
surrounding envelope, which acts as an oven (see van der Tak et al. 1999). Assuming a temperature 
of 200 K for the compact dust, a lower limit to the mass of ~ 10 Mq is obtained. In the case of 
NGC 7538 IRSl, where free-free and dust emission both contribute at these frequencies (Woody 
et al. 1989; Akabane et al. 1992), the physical temperature may be > 100 K, and the beam dilution 
correspondingly more. For the other sources, the compact emission is most likely due to dust, 
perhaps in the form of a dense shell or of a circumstellar disk. The implied radius of the emitting 
region, < 300 — 600 AU at 2 — 4 kpc, seems small to hide an entire outflow lobe, but presumably, 
the dense part of the flow which emits in high- J CS lines is confined to its center. 



6.4. Comparison with Low-mass Objects 

This section discusses possible origins of the spread in a within our sample. First, the value 
of a may be related to other physical parameters, in particular luminosity or envelope mass. The 
structure of the material surrounding pre-main sequence stars of low ('^ solar) and intermediate 
(~ 3 — 8 solar) mass has received considerable attention in the recent literature: e.g. Ladd et al. 
(1991); Butner et al. (1991); Natta et al. (1993); Butner et al. (1994); Hogerheijde et al. (1999). 
These studies generally find a 2 it 0.3 inside radii of 0.1 pc, masses of <10 Mq and mean densities 
of 10^ cm""^. Of these properties, only the radii are similar to those of the objects studied here; 
the masses and mean densities arc at least two orders of magnitude smaller. The density gradients 
are significantly steeper than found in this paper, but the gradients for intermediate-mass stars 
arc also lower than those for low-mass stars (Natta et al. 1993). Indeed, in an earlier study of 
regions of high-mass star formation, Zhou et al. (1994) suggested that more massive star-forming 
regions tend to have flatter density distributions. However, Figs. 13a and b show that within our 
sample, a does not appear correlated with source luminosity or envelope mass. We did not find a 
relation with any other physical parameters either, such as turbulent pressure as measured by the 
line width. 

Power laws have been proposed for the density structures of the envelopes of young stars, with 
the index depending on the dominant term in the pressure. If support against gravitational collapse 
is primarily due to thermal pressure, a value of a = 2 is expected, while if the dominant pressure is 
of nonthermal origin, the density structure should follow an a = 1 law. In between these two cases, 
a continuum of solutions can be constructed (Lizano & Shu 1989; Myers & Fuller 1992; McLaughlin 
&; Pudritz 1997). Our observations may thus indicate that the cores where massive stars form differ 
from those where low-mass stars form in that they are supported against collapse by a different 
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mechanism. Within this interpretation, the spread in a may be due to source evolution from an 
(a = 1.0) logatrope to a collapsing region (a = 1.5). 

6.5. Tracing Envelope Evolution 

Despite having similar radio and infrared properties oiir sources show different degrees of 
dispersal and warming-up of their envelopes. This is shown by Boogert et al. (2000) and van 
Dishoeck (1998) for a source sample similar to ours based on several tracers: the far-infrared color, 
the fraction of crystalline CO2 ice, the gas/solid ratios of CO, CO2 and H2O, and the excitation 
temperature of CO. These quantities all trace the envelope temperature, but all in a different way: 
some trace the dust, others the gas, and some trace small scales and others large scales. As shown 
by Boogert et al., the temperature indicators correlate with each other. Of particular interest is 
the fraction of crystalline CO2 ice, because crystallization is an irreversible process, so that this 
indicator measures the maximum temperature, while the others measure the current temperature. 
Hence, the temperature variations are not random fluctuations (like in FU Orionis objects), but 
reflect a systematic increase with time of the temperature of the envelope as a whole. 

Figure 14a plots one of these indicators, the far-infrared color, versus the ratio of M(< tq) 
to stellar mass, measured as L^/^-^. This quantity is equivalent to the ratio of submillimeter to 
bolometric luminosity, often used to trace the evolution of young low-mass objects (Andre et al. 
1999). The color i^45/i^ioo is the ratio of the flux densities at 45 and 100 /im, both observed with 
ISO-LWS in an 80" beam, which is large enough to cover the entire envelope. These data were 
provided by A. Boogert (1999, priv. comm.). The plotted quantities are seen to be correlated, with 
higher temperatures (traced by a larger F45/F100) corresponding to lower values of M(< ro)/L^f^'^. 
Thus, the temperature variations are indeed due to dispersal of the circumstellar material and can 
be used to trace evolution. 

Interestingly, a is not seen to be correlated with far-infrared color (Fig. 14b) , and we conclude 
that on the time scale of the dispersal of the envelopes of massive young stars, the density structure 

of the envelopes does not change significantly. The masses of the compact sources detected by us 
and by Woody et al. (1989) and Wyrowski et al. (1999) do not show a relation with M(< ro)/L^^^'^ 
either, although the uncertainty in these masses is rather high. 

In addition to M(< ro)/iv^/'^'^, we have considered two "derivative" evolutionary indicators. 
The first is the bolometric temperature Tboi, introduced for low- mass objects by Myers & Ladd 
(1993). As the envelopes in those systems are dispersed (probably by the bipolar outflows), their 
bolometric temperatures rise because the far-infrared emission decreases while the near-infrared 
emission increases (Myers et al. 1998). This process makes Tboi increase monotonically from ^ 30 K 
to 3000 K on a time scale of 10^ yr. For all but two of our objects, values of Tboi based on the 
spectral energy distributions of Fig. 6, are 50 — 150 K. This small spread suggests that all the 
sources are at similar evolutionary stages; Tboi may also work less well than F45 / Fioo because near- 
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and mid-infrared emission have a dependence on source orientation. 

Another potential evolutionary indicator is the radio continuum flux, which for a given lumi- 
nosity and distance is expected to increase with time as the dusty envelope is cleared away and the 
Lyman continuum photons can all go into ionizing the gas without being absorbed by dust. Wc 
have normalized the radio continuum flux densities in Tabic 1 to the value expected for an H II 
region ionized by a main sequence star of the same luminosity and at the same distance as the 
source, following Churchwell (1993). However, Figs. 14c and d show that neither the radio con- 
tinuum emission nor the bolometric temperature correlates well with M(< ro)/L^^^'^, suggesting 
these are less useful evolutionary indicators for very young massive objects than Kis/Fioo- Possibly 
the radio emission arises in a wind (§ 2), or the stellar surface temperatures are signiflcantly below 
main sequence values. 

We conclude that dispersal of the envelopes underlies the evolution of both high-mass and 
low-mass objects, but that the observational appearance is very different. The dispersal process 
can be understood in more detail by recalling that the brightest mid-infrared sources in our sample 
are the weakest at radio wavelengths and vice versa (Table 1). This anticorrelation suggests that 
the erosion of the envelope proceeds from the inside out. Massive stars are efflcient at removing 
circumstellar dust by their ionizing UV radiation and by their strong winds. When the innermost, 
hottest region of dust disappears, the mid-infrared emission will decrease, but the radio continuum 
will increase because a larger region can be ionized. The compact dust component detected with 
OVRO is optically thick in the bright mid-infrared sources but thin in the others, which trend 
may also trace dispersal of hot dust. The far-infrared emission will not be affected since it arises 
on larger scales. These trends are consistent with the properties of our sample, although other 
processes may play a role as well. 



7. Summary and Conclusions 

We have presented maps and spectra of 14 regions of massive star formation in continuum and 
molecular lines at submillimeter wavelengths. The data are used to develop models of the physical 
structure of the circumstellar envelopes on 10^ — 10^ AU scales. The column density is derived 
from the submillimeter continuum flux densities, and the temperature structure is calculated self- 
consistently using the size of the CS emission and the sources' luminosities and distances from the 
literature. The density structure is constrained with emission lines of CS, C^^S and H2CO. The 
following main conclusions are found: 

1. The physical structure of the envelopes of deeply embedded massive young stars is charac- 
terized by radii of 3 — 910^ AU, masses of 10^ — 10^ Mq inside these radii and visual extinctions of 
~ 10^ — 10"^ magnitudes. Temperatures increase from ~ 20 K at the outer edge to several 100 K 
at 10^ AU from the star; densities increase from ~ 10^ cm~^ to ~ 10^ cm~^. The slope of 
the density gradient a is 1.0 — 1.5, signiflcantly smaller than the value of « 2 usually found for 
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low-mass objects. This difference in density structure may be due to nonthermal pressure resisting 

gravitational collapse, while thermal pressure dominates in lower-mass objects. An a = 2 density 
structure is found here for two hot cores, but this is an upper limit since the CS abundance may 
be enhanced in the warm gas close to the star. 

2. The shapes of the envelopes deviate somewhat from spherical, as shown by modeling infrared 
absorption line data of ^^CO, the near-infrared continuum and the CS line profiles. The data are 
consistent with a decrease of the optical depth by a factor of pa 3 in the central $ 10" area. The 
brightness of our main sample at mid-infrared wavelengths is partly due to this optical depth effect, 
but also to the flatter density structure. Inhomogeneities on top of the power law structure are 
small, since the masses obtained by integrating the power law models agree with masses found 
from the virial theorem to a factor of 3. The values of a found from CS are verified by model fits 
to the maps of dust emission which trace the column density distribution. The only exception is 
S 140 IRSl, probably because our assumption of a single central heating source is invalid in this 
case. 

3. Modeling of C^'^0 emission lines shows that pa 40 — 90% of the CO gas is depleted onto 
dust, assuming a dust opacity of pa 1 cm~^ g~^ at A1300/xm, as found previously for one of our 
sources, GL 2591 (van der Tak et al. 1999). The derived CO abundances correlate well with the 

mass-weighted temperature in our models. This result suggests that in these sources, freeze-out 
and thermal dcsorption control the gas-phase abundance of CO. The CS abundance is 3 x 10~^ on 
average, ranging from 0.4 x 10~^ in the cold source GL 7009S to 1 — 2 x 10~^ in the "hot core" -type 
objects. 

4. Dense outflowing gas is seen in the CS and H2CO line wings, but much more so at blue- than 
at redshifted velocities. This asymmetry may be due to absorption by dense, optically thick material 
within 10" from the star. The opacity of the power law models is a factor of 10 — 100 below that 
required for this absorption. However, interferometric continuum observations at A1300 — 3500)um 
show compact emission, probably from an 0^'3 diameter optically thick dust component. The 
emission is a factor of 10 — 100 stronger than expected for the envelopes seen in the single-dish 
data, and may be due to a dense shell or a circumstellar disk. This component may explain the 
asymmetric CS and H2CO line profiles. 

5. The evolution of these high-mass sources is traced by the overall temperature (measured 
by the far- infrared color), which increases systematically with decreasing ratio of envelope mass to 
stellar mass. The observed anticorrelation of near-infrared and radio continuum emission suggests 
that the disruption of the envelope proceeds from the inside out. Conventional tracers of the 
evolution of low-mass objects do not change much over this narrow age range. We conclude that 
the cvohition of high-mass and low-mass envelopes have the same underlying mechanism (envelope 
dispersal), despite their different observational properties. 
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Fig. 1. — Line profiles of CS and C^^S, observed with the JCMT and IRAM 30m telescopes. For 
clarity, the IRAM data have been divided by 2 and the C'^'^S data have been multiplied by 2 for 
the hot cores and by 3 for the other sources. At blueshifted velocities, wings are visible on the CS 
lines, most prominently in the IRAM data. 
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Fig. 2. — Column densities of CO derived from C^''0 emission compared to values derived from 
infrared absorption of ^^CO. The dotted line indicates the relation expected for a constant-density 
model. 
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Fig. 3.— CSO maps of the C^^O J = 2 ^ 1 and CS J = 5 ^ 4 and 7 ^ 6 lines. Contours are 10% 

of the peak flux, starting at 10% for CS and 30% for C^^O. Integration intervals and peak C^^O 
fluxes are: 19 to 26 km s^^ and 3.5 K km s^^ for GL 2136, -15 to -10 km s^^ and 4.0 K km s"^^ for 
GL 490, -18 to km s^^ and 7.7 K km s^^ for NGC 6334, -65 to -50 km s'^ and 5.0 K km s'^ for 
NGC 7538, -12 to -2 km s^^ and 3.5 K km s'^ for S 140, 30 to 44 km s"^ and 3.8 K km s"^ for 
W 33A, -58 to -38 km s'^ and 5.1 K km s"^ for W 3 (H2O) and -50 to -32 km s"^ and 5.2 K km s"^ 
for W 3 (Main). See Table 2 for the peak fluxes of the CS lines. 
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850 iim, observed with SCUBA on the JCMT. Contours are 10 to 90% of the peak brightness, 
in steps of 10%. Peak brightness (in Jy/beam) is 45.8 for W 33A, 7.7 for W 3 IRS5, 27.7 for 
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Fig. 5.— OVRO continuum maps of GL 2136, W 33A, NGC 7538 IRSl and NGC 7538 IRS9. 
Contour levels (in mJy/beam) are: 6, 18, 30, 42 (a); 8, 16, 24, 32, 40, 48, 56 (b), 40, 80, 120, 160, 
200, 240, 280 (c); 6, 18, 30, 42 (d), 50, 550, 1050, 1550 (e), 500, 2000, 3500, 5000, 6500, 8000 (f). 
The beam sizes and shapes are drawn in the bottom right corner of each panel. 
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Wavelength (/U.m) 



Fig. 6.— 

Observed continuum spectra (open circles) and models for several values of a. The symbol legend 
is in the top left panel. For W 33A, a model using grains without ice mantles is also shown. Data 
are taken from Campbell et al. (1995) for W 3 IRS5; Chini et al. (1991) for GL 490; Moorwood & 
Salinari (1981), Evans et al. (1979), Dyck & Simon (1977), Stier et al. (1984), Jaffe et al. (1984) 
and Cheung et al. (1980) for W 33A; Allen et al. (1977), Lebofsky et al. (1976), van Dishoeck 
(priv. comm.) and Kastner ct al. (1994) for GL 2136; Dartois et al. (1998a) and McCutcheon et al. 
(1995) for GL 7009S; Evans et al. (1989), Willner et al. (1982) and Zhou et al. (1994) for S 140 IRSl, 
Werner et al. (1979) for NGC 7538 IRSl; Werner et al. (1979) and G. Sandeh (priv. comm.) for 
NGC 7538 IRS9; Wilking et al. (1990) and Cesaroni et al. (1997, 1999) for IRAS 20126+4104; 
Thum &: Lemke (1975), Harvey et al. (1986, 1977), Richardson et al. (1989) and Chandler et al. 
(1993) for DR 21 (OH); Wynn-Wilhams et al. (1972), Keto et al. (1992) and Thronson & Harper 
(1979) for W 3 (HgO), Straw & Hyland (1989) and Harvey & Gatlcy (1983) for NGC 6334 IRSl; 
and Moorwood & Salinari (1981), Lightfoot et al. (1984) and Harvey et al. (1994) for W 28 A2. 
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Figure 6 (continued) 



-40- 




Fig. 7. — Fit quality parameter ^ plotted as a function of density law exponent a and CS abun- 
dance. Contours increase by 1 and start at 1 for W 3 IRS5, GL 7009S, IRAS 20126, W 3 (H2O) 
and NGC 6334 IRSl, 3 for NGC 7538 IRSl and NGC 7538 IRS9, and at 2 for the other seven 
sources. The best-fitting model is indicated by a star. In the top right corner, the number of CS 
and C^^S lines included in the calculation is listed. 
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Figure 7 (continued) 
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Fig. 8. — Radial profiles of 350 /im emission observed with SHARC (a) and SCUBA (b). Dotted 
lines are model fits from the power law models for a = 0.5, 1.0, 1.5 and 2.0 (top to bottom). The 
solid line is the model that fits the CS excitation best. 
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Figure 8 (continued) 
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Fig. 9. — Observed and modeled column density in the rotational states of ^^CO up to J = 25 for 
the NIR-bright sources. Data are from Mitchell et al. (1990) except GL 7009S where they are from 
Dartois et al. (1998). The model values have been scaled to the observed total column density. 
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Fig. 10. — Comparison of cloud masses derived from the power law models with those derived from 
the virial theorem. The quantity Upc is the density at a distance of 1 pc from the center; AV is the 
line width from Table 3 and the other quantities are in Table 5. 
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Fig. 11. — Abundance of CO derived from ^CO observations plotted against the mass- weighted 
temperature in our models. 
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Fig. 12.— 



Continuum visibilities observed with OVRO compared to the power law models. 
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Fig. 13. — Slope of density gradient a versus source luminosity and mass. Filled symbols indicate 
the main sample, open symbols the comparison sample. The plotted quantities are explained in 
the text. 
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Fig. 14. — Possible evolutionary indicators plotted versus ratio of envelope mass to stellar mass: 
far- infrared color (a), slope of density gradient a (b), normalized radio continuum flux density 
(c) and bolometric temperature (d). Filled symbols indicate the main sample, open symbols the 
comparison sample. See the text for details. 
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Table 1. Source sample. 



IRAS PSC 


Name 


RA (1950) 


Dec (1950) 


Luminosity 


Distance 


Si, (6 cm) 


References 
















(L©) 


(kpc) 


(mJy) 












MIR-bright massive 


YSOs 








02219+6152 


W 3 IRS5 


02'' 


21m 


53? 1 


+61° 52' 


20" 


1.7 X 10^ 


2.2 


10.7 


1,2,21 


03236+5836 


GL 490 


03'' 


23"" 


38! 9 


+58° 36' 


33" 


2.2 X 10^ 


1 


3.2 


3,4,22 


18117-1753 


W 33A 


18'' 


11™ 


43=7 


-17° 53' 


02" 


1 X 10"' 


4 


1.9 


5,23 


18196-1331 


GL 2136 


IS'' 


19"' 


36? 6 


-13° 31' 


40" 


7 X 10'' 


2 




6 


18316-0602 


GL 7009S 


18'' 


31™ 


41?6 


-06° 02' 


35" 


3 X 10"' 


3 


2.7 


7 


20275+4001 


GL 2591 


20'' 


27m 


35? 8 


+40° 01' 


14" 


2 X 10" 


1 


0.4 


8,9,24 


22176+6303 


S 140 IRS 1 


22'' 




41?08 


+63° 03' 


41'.' 6 


2 X 10" 


0.9 


5.8 


10,11,25 


23116+6111 


NGC 7538 IRS 1 


23'' 


lim 


36?7 


+61° 11' 


50'.'8 


1.3 X 10^ 


2.8 


111 


12,13,26 


23118+6110 


NGC 7538 IRS 9 


23'' 




52? 8 


+61° 10' 


59" 


4.0 X 10* 


2.8 


< 0.5 


12,13,23 










MIR-wcak massive 


YSOs 








20126+4104 


IRAS 20126+4104 


20'' 


12™ 


41.0^ 


+41° 04' 


21" 


4.5 X 10^ 


1 


< 0.3 


14,9,27 




DR 21 (OH) 


20'' 


37™ 


14.2'' 


+42° 12' 


11" 


1 X 10^ 


1 


< 10 


15,9,28 


02232+6138 


W 3 (H2O) 


02'' 


23™ 


17.3^ 


+61° 38' 


58" 


2 X 10* 


2.2 


1.5 


16,2,29 


17175-3544 


NGC 6334 IRS 1 


17'' 


17™ 


32.0* 


-35° 44' 


05" 


1.1 X 105 


1.7 


1200 


17,18,30 


17574-2403 


W 28 A2 


17'' 


57™ 


26.8* 


-24° 03' 


57" 


1.3 X 10^ 


2.0 


2100 


19,20,31 




(=G5.89-0.39) 





















Note. — The first reference is for the luminosity, the second for the distance, and the third for the radio data. Spectrophotometric 
distances are given to 1 decimal place, and the corresponding luminosities are accurate to a factor of 2. The other distances are 
kinematic, in which case the luminosity is uncertain by a factor of 4. The kinematic distances to W 33A and GL 2136 were 
derived from data presented in this paper. 

References. — (1) Ladd et al. 1993; (2) Humphreys 1978; (3) Chini et al. 1991; (4) Harvey et al. 1979, but see Snell et al. 1984; 
(5) Giirtler et al. 1991; (6) Kastner et al. 1994; (7) McCutcheon et al. 1995; (8) Lada et al. 1984; (9) Distance discussed in van 
der Tak et al. 1999; (10) Lester et al. 1986; (11) Crampton & Fisher 1974; (12) Werner et al. 1979; (13) Crampton et al. 1978; 
(14) Cesaroni et al. 1997; (15) Chandler et al. 1993; (16) Turner & Welch 1984; (17) Harvey & Gatley 1983; (18) Neckel 1978; 
(19) Harvey et al. 1994; (20) Acord et al. 1998; (21) Tieftrunk et al. 1997; (22) Simon et al. 1983; (23) Rengarajan & Ho 1996; 
(24) Campbell 1984; (25) Evans et al. 1989; (26) Pratap et al. 1992; (27) Tofani et al. 1995; (28) Johnston et al. 1984; (29) Reid 
et al. 1995; (30) Rodriguez et al. 1982; (31) Wood & Churchwell 1989 
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Table 2. Observed line fluxes / T^\,dV (K km s-^). 



Line 


W 3 


GL 490 


W 33A 


GL 2136 


S 140 NGC 7538 


NGC 75.38 


NGC 6334 


W 3 




IKSr, 








IK SI 


IHSl 


1KS9 


IRSl 


(il,0) 










JCMT observations 










Cl^O 2-l('') 


12.2 




11.3 


6.8 


6.4 


8.7 


3.8 


26.5 




3-2 


16.6 


4.9 


15.7 


9.9 


10.8 


17.9 


4.5 


57.9 


15.1 


CS 5-4 




16.5 


36.2 


11.2 


30.2 


54.8 


20.6 






7-6 


25.2 


8.0 


17.4 


4.7 


22.3 


56.0 


20.5 


162. 


83.7 


10-9 


16.1 










34.4 


3.3 




78.4 


C^^S 5-4('') 


3.7 


1.1 


5.6 


1.1 


2.7 


7.0 


2.4 


56.5 


20.7 


7-6 


2.5 


0.6 


5.3(rf) 


< 0.4 


2.1 


7.9('*) 


3.1 


64.8('*) 


25.8 


10-9 


2.2 










4.6 


< 3.4 




31.6 


rdoCjO CiiYi - — zno 


8.5 


4.4 


13.1 


4.4 


13.7 


27.0 


14.2 


63.2 


20.2 


3i2 — 2ii 


19.0 




15.0 


6.3 


19.2 


97 n 


13.6 


63.0 


41.2 


322 — 221 


3.6 


0.7 


3.5 


1.6 


3.5 


7.1 


3.7 


26.2 


6.9 


505 - 4o4 


5.2 




10.9 


4.3 


8.7 


18.2 


9.5 


82.9 


31.1 


5l5 — 4l4 


14 




20.0 


7.4 


20.5 


28.8 




1 




524 — 423 


2.7 


1.6 


4.5 


2.2 


3.7 


7.9 


4.5 


52.9 


19.5 


532 - 431 


2.9 


1.5 


5.4 


1.5 


3.6 


8.8 


4.5 


48.1 


33.4(«: 


533 - 432 


2.6 


1.2 


4.9 


2.6 


3.3 


9.9 


4.2 


50.6 


24.7 


542/41 — 441/40 


0.9 


< 0.7 


1.5 


< 0.4 


< 0.5 


2.4 


1.3 


29.4 


8.3 


7l7 - 6l6 


21.9 








7.3 


15.8 


10.8 




24.3 








IRAM 30m observations 










CS 5-4 






80.7 


22.5 


37.7 


79.0 


37.2 






C'^-^S 2-1 






7.6 


1.8 


3.6 




4.1 






C^^S 3-2 








2.4 


4.3 




5.4 














CSO observations 










CS 5-4('=) 


21.5 




20.3 


8.3 


30.5 


54.4 


18.2 


78.6 


33.5 


CS 7-6 


20.0 


4.9 


10.5 


6.6 


20.7 


50.7 


< 0.2 


103.6 


64.2 


C^^S 7-6 


















15.8 



("^'For this line, a flux of 2.9 K km s"! was measured for IRAS 20126 and 24.5 K km s^^ for DR 21 (OH). 

('')For this hne, a flux of 0.5 K km s"! was measured for GL 7009S, 3.2 K km s"! for IRAS 20126 and 57.9 K km s"! for 

W 28 A2. 

('=)For this line, a flux of < 0.8 K km s"! was measured for GL 7009S, < 0.6 K km s'^ for IRAS 20126, 59.7 K km s'^ for 
DR 21 (OH) and 123.2 K km s"! for W 28 A2. 

('^^Blend with triplet of CN lines; fit constrained by requiring same Ay as other C**S and '^^CO lines. 

('^)Blend with CH3OCHO line complex; fit constrained by requiring same Ay as other H2CO lines. 
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Table 3. Single-dish line emission: basic parameters. 



Source 


Vlsr 


AV 


N(CO) 


D(CS) 




km s"-^ 


km 


s-i 


10^9 cm-2 


arcsec 


W 3 IRS5 


-38.4 (0.3) 


2.7 


(1.1) 


3.7 


55 


GL 490 


-13.3 (0.2) 


3.3 


(0.6) 


0.8 


38 


W 33A 


-1-37.5 (0.9) 


5.4 


(0.5) 


3.1 


37 


GL 2136 


+22.8 (0.1) 


3.1 


(0.4) 


1.4 


36 


GL 7009S 


+40.3 (0.6) 


5.0 


(1.0) 






GL 259l(") 


-5.5 (0.2) 


3.3 


(0.6) 


3.4 


52 


S 140 IRSl 


-7.0 (0.2) 


3.3 


(0.3) 


1.4 


64 


NGC 7538 IRSl 


-57.4 (0.5) 


4.1 


(1.4) 


3.9 


52 


NGC 7538 IRS9 


-57.2 (0.3) 


4.1 


(0.5) 


1.0 


47 


IRAS 20126('') 


-3.7 (0.2) 


3.2 


(0.2) 


0.8 


40 


DR 21 (OH) 


-3.1 (0.2) 


4.5 


(0.2) 


15.6 


56 


W 3 (H2O) 


-47.6 (0.6) 


5.8 


(0.6) 


8.5 


41 


NGC 6334 IRSl 


-7.4 (0.2) 


5.3 


(0.3) 


16.9 


50 


W 28 A2f'') 


+10 




5.7 




38 



Note. — The CO column densities in column 4 arc derived from 
C^^O emission assuming ^^0/^^0=2500, refer to a 14" beam and 
have a 30% calibration uncertainty. Column 5 lists the half-power 
diameters of the CS J = 5 ^ 4 emission mapped with the CSO. 

("^Size of CS emission from data presented in Carr et al. (1995) 

(^)Size of CS emission from data presented in Cesaroni et al. 
(1997) 

("^^Velocity and line width taken from Plume et al. (1997) 
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Table 4. Results of OVRO observations. 



Source 


Frequency 
(GHz) 


RA (1950) 


Dec (1950) 


Flux Density 
(mJy) 


Beam Size 
(arcsec) 


GL 2136 


86 


18:19:36.62 (0.02) 


-13:31:43.8 (0.3) 


61 


3.0 X 2.8 


W 33A MMl 


86 


18:11:44.22 (0.02) 


-17:52:58.2 (0.2) 


24 


3.3 X 2.7 




106 






54 


2.9 X 2.6 




233 






190 


2.0 X 1.6 


W 33A MM2 


86 


18:11:43.91 (0.02) 


-17:52:59.8 (0.2) 


38 


3.3 X 2.7 




106 






30 


2.9 X 2.6 




233 






140 


2.0 X 1.6 


NGC 7538 IRSl 


86 


23:11:36.64 (0.02) 


+61:11:49.8 (0.2) 


1500 


2.6 X 2.0 




107 






1800 


2.4 X 1.9 




115 






2000 


2.5 X 1.7 




230 






3500 


1.1 X 0.9 


NGC 7538 IRS9 


107 


23:11:52.90 (0.07) 


+61:10:59.2 (0.5) 


43 


2.4 X 2.0 




115 






95 


2.4 X 1.9 



Note. — 



The flux densities have a calibration uncertainty of « 30%. 
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Table 5. Parameters of Best Fit Models 



Source 


^0 


a 


no 


M(< ro) 


My 


T 


CO/H2 


CS/H2 


H2CO/H2 




10=^ AU 




10^ cin-=^ 


Mo 


Mo 


K 


xlO^ 


xlO*^ 


xlO'' 


W 3 IRS5 


60 


1.5 


2.3 


262 


355 


33 


1.6 


5. 


3. 


GL 490 


19 


1.5 


8.6 


30 


168 


19 


0.4 


1. 


1. 


W 33A 


74 


1.0 


6.8 


1089 


1984 


20 


0.5 


5. 


4. 


GL 2136 


35 


1.25 


3.6 


72 


316 


28 


1.2 


4. 


8. 


GL 7009S 


91 


0.5 


17. 


3955 


2218 


21 




0.4-0.8 




GL 2591 


27 


1.0 


5.8 


44 


268 


28 


1.5 


10. 


4. 


S 140 IRSl 


29 


1.5 


3.6 


46 


128 


27 


1.0 


5. 


5. 


NGC 7538 IRSl 


72 


1.0 


5.3 


815 


1128 


25 


0.6 


10. 


10. 


NGC 7538 IRS9 


66 


1.0 


3.9 


430 


1016 


20 


0.3 


10. 


10. 


IRAS 20126 


39 


1.75 


2.2 


84 


294 


27 


0.5 


3. 




DR 21 (OH) 


29 


1.75 


23. 


350 


429 


15 


0.8 


3. 


5. 


W 3 (H2O) 


45 


2.0 


5.3 


385 


932 


26 


1.0 


10. 


3. 


NGC 6334 IRSl 


43 


2.0 


6.2 


382 


725 


35 


1.5 


12. 


7. 


W 28 A2 


37 


1.5 


9.9 


292 


982 


42 


1.7 


20. 





Note. — The parameters in columns 2-4 specify the density structure in the form n(r) = no(r/ro)~", 
with the reference radius ro half the diameter listed in Table 3. Column 5 gives the integrated mass enclosed 
within ro, and column 5 the virial mass within ro- 



